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1.  Introduction 


Conventional  methods  for  discriminating  between  earthquakes  and  explosions  at  re¬ 
gional  distance  have  concentrated  on  extracting  specific  features  from  the  waveforms  of 
the  P(usually  Pg)  zmd  S  (usually  Lg)  phases.  The  specific  features  considered  generally 
are  zunplitude  ratios,  measures  of  waveform  complexity  or  various  kinds  of  spectral  ratios, 
suggesting  that  the  main  characterization  of  the  differences  between  earthquakes  and  ex¬ 
plosions  reduces  to  differences  between  the  spectra  or  differences  between  the  waveforms. 
0\ir  objective  here  is  to  compare  some  of  the  classical  discriminants  in  the  literature  with 
two  methods  based  on  statistical  optimality  criteria.  We  consider  first  an  optimal  nonpara- 
metric  discriminator,  derived  under  the  assumption  that  the  earthquake  and  explosion  P 
and  S  phases  are  uncorrelated  stationary  Gaussian  processes  with  unequad  spectra.  The 
likelihood  criterion  that  obtains  displays  the  optimal  statistic  as  the  result  of  comparing 
spectral  matches  between  the  observed  series  and  the  average  earthquake  spectrum  against 
a  comparable  match  with  the  explosion  sp>ectrum.  Two  vauiations  based  on  information 
theoretic  principles  are  also  investigated.  A  parametric  discriminator  models  the  P  phaise 
as  a  modulated  autoregressive  process,  where  the  modulating  function  is  consistent  with 
models  for  earthquake  and  explosion  waveforms  foimd  in  the  literature.  The  nonparamet- 
ric  method  is  obviotisly  tuned  to  spectral  differences  whereas  the  parametric  method  is 
closely  adlied  with  notions  relating  to  complexity  and  amplitude  ratios. 

Niunerous  investigators  have  pointed  out  that  the  logarithms  of  Pg/Lg  amplitude 
ratios  tend  to  be  lower  for  earthquakes  than  for  explosions  (see,  for  example  Blandford, 
1981,  Bennett  and  Murphy,  1986,  Taylor  et  al,  1989).  This  idea  has  been  extended  to 
include  a  consideration  of  spectral  ratios  involving  the  P  and  S  groups.  Bennett  and 
Murphy  (1986)  note  that  for  western  U.S.  events,  earthqualce  Lg  spectra  contained  more 
high  frequencies,  and  that  the  ratio  of  the  logarithms  of  low  frequency  (.5-1  Hz)  Lg  to 
higher  frequency  Lg  (2-4  Hz)  tend  to  be  larger  for  explosions.  Taylor  et  al  (1989)  also  use 
this  ratio  over  the  frequency  bands  (1-2  Hz)  and  (6-8  Hz)  and  extended  the  consideration 
to  the  Pg  phase.  Dysart  and  PuUi  (1990, 1992)  have  edso  considered  various  spectral  ratios 
P/S  for  Scandinavian  events  and  have  developed  neural  networks  as  an  alternative  to 
simple  linear  combinations  of  features  for  discrimination.  They  note  that  the  P/ S  spectral 
ratios  are  generally  higher  for  explosions  than  for  earthquakes.  Finally,  Kim  et  al  (1992) 
note  that  for  eastern  U.S.  events  the  ratios  of  Pg  to  Lg  spectra  are  generally  higher  for 
explosions.  Some  early  results  using  the  three  coefficients  in  an  third-order  autoregressive 
model  for  the  coda  as  features  are  available  in  Tjpstheim  (1974). 
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The  spectral  methods  discussed  above  generate  features  that  can  differ  for  earthquakes 
and  explosions  in  certain  specific  data  bases.  Such  features  as  P  to  S  amplitude  ratios 
or  spectral  ratios  within  phases  or  between  phases  can  then  be  put  into  a  vector  and 
transformed  (usually,  logarithms  are  taken  for  a  better  approximation  to  normality)  in 
order  to  apply  one  of  the  standard  linear  or  quadratic  discriminant  analysis  techniques. 
Examples  are  Shumway  and  Blandford  (1974),  Taylor  et  al  (1989),  Dysart  and  Pulli  (1990), 
PuUi  and  Dysart  (1992)  and  Kim  et  al  (1992).  Two  features  at  a  time  are  often  plotted  in 
various  combinations  for  earthquakes  and  explosions  to  show  graphically  the  sepauration  of 
the  two  populations  (see  the  above  references  and  also  Bennett  and  Murphy,  1986). 

Shumway  and  Blandford  (1974)  introduced  am  optimal  method  combining  optimal 
linear  and  quadratic  discriminant  functions.  The  criterion  was  based  on  modeling  the 
imderlying  short  period  teleseismic  P  waveforms  as  Gaussiam  processes  differing  in  both 
the  mean  vailue  signals  and  the  spectrad  densities.  For  regional  events,  it  is  clear  that  the 
notion  of  a  fixed  mean  P  or  S  waveform  which  differs  for  earthquaJces  and  explosions  is 
not  a  relevant  comparison  but  that  the  notion  that  the  eau-thquaJces  and  explosions  differ 
only  in  their  spectra  does  make  a  lot  of  sense.  This  implies  that  the  quadratic  part  of  the 
optimal  detector  used  by  Shumway  and  Blandford  (1974)  (see  also  Shumway,  1982,  1988) 
will  have  the  lowest  misclassification  rate  of  any  function  based  on  the  spectra,  including 
those  given  in  the  preceding  paragraph.  Such  a  detector  has  been  applied  several  times 
in  the  literature  to  seismic  discrimination,  by  other  investigators  as  well  (see,  for  example 
Deirgsdii-Noubary  and  Laycock,  1981,  Alagon,  1989).  We  develop  amd  apply  a  modified 
version  of  such  a  detector  to  P  and  S  phases  from  population  of  regional  Scamdinavian 
earthquakes  and  explosions  given  in  Blandford  (1993).  Since  the  model  depends  only  on 
stationarity  of  the  series  and  not  on  a  specific  peurametric  model  for  the  spectrum,  we  refer 
to  it  as  the  optimal  nonparametric  discriminator. 

Modifications  to  the  nonparametric  version  of  the  likelihood  detector  can  be  made 
baised  on  information  theoretic  principles.  For  example,  the  optimum  quadratic  detector 
has  am  information  theoretic  interpretation  in  terms  of  the  minimum  discrimination  infor¬ 
mation  statistic  (MDI),  of  KuUback  (1978).  Such  a  discriminant,  defined  as  the  difference 
of  the  discrepancies  between  the  sample  spectnim  of  the  event  to  be  classified  and  the  the¬ 
oretical  earthquake  and  explosion  spectra,  has  excellent  theoretical  properties  as  discussed 
by  Zhang  and  Taniguchi  (1992,  1993).  They  suggest  am  alternate  discriminamt  that  is  ro¬ 
bust  to  peak  contamination  which  is  based  on  the  Renyi  index  of  index  a  as  discussed  by 
Parzen  (1990)  (see  Renyi,  1961).  Zhan  and  Taniguchi  call  the  discriminant  the  a-entropy. 
We  will  apply  the  MDI  and  Zhang- Taniguchi  (ZT)  discriminauits  to  the  earthquake  and 
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explosion  popvdations  and  show  that  the  ZT  discriminant  offers  several  2ulvantages  over 
the  conventional  likehhood  discriminant. 


Blandford  (1993)  discusses  a  notion  of  complexity  as  a  discriminant  and  shows  that 
it  has  potential  for  discrimination  of  the  events  in  the  Scandinavian  database.  The  idea 
relates  to  the  often  quoted  statement  by  analysts  that  complexity  is  a  strong  component 
of  their  visual  procedure  for  discrimination.  Blandford  proposes  a  notion  of  complexity 
related  to  the  observation  that  explosions  generate  usually  am  impulsive  signal  whereas 
earthquakes  tend  to  generate  a  more  emergent  signal.  We  develop  here  a  parametric 
discriminator  by  assuming  that  the  earthquake  and  explosion  populations  can  be  expressed 
as  uniformly  modulated  autoregressive  process;  parauneters  of  the  modulating  function 
characterize  separately  the  emergent  and  impulsive  properties  of  the  signal.  Our  \mderlying 
model  for  complexity  is  taken  directly  from  a  suggestion  of  Daurgahi-Noubary  (1992)  that 
is  based  on  standard  source  theory. 


2.  Data  Compilation 

For  our  test  data,  we  use  a  subset  of  stations  recording  8  earthquakes  and  8  explosions 
in  Scamdinavia  from  the  arrays  NORESS,  ARCESS  and  FINESS  as  described  in  Blandford 
(1993).  According  to  Blandford,  “The  events  were  selected  with  consideration  for  having 
sxiflBcient  S/N  at  single  elements  so  that  all  phases  could  be  clearly  seen  on  ail  components 
of  a  single  instrument  •  •  •”.  From  Table  1  in  Blandford,  we  took  explosions  three  through 
ten  and  from  Table  2,  we  used  all  eight  earthquakes.  All  events  chosen  by  Blandford 
were  on  or  near  land  and  were  distributed  xmiformly  over  Scandinavia  to  minimize  the 
possibility  that  discriminators  might  be  keying  on  location  or  land>sea  differences.  Figure 
1  shows  portions  of  typical  earthquakes  zind  explosions  (sampled  at  20  Hz)  along  with  the 
portions  of  the  record  that  we  visually  determined  as  the  P  and  S  phases.  We  did  not 
identify  specific  phases  through  velocity  computations  but  simply  chose  fairly  broad  (25 
second)  windows  that  seemed  to  include  the  P  and  S  phases. 

Qualitatively,  we  note  that  the  earthquake  heis  a  much  smaller  amplitude  P/S  ratio 
than  the  explosion  and  we  note  the  relatively  complex  P  phase  which  tends  to  be  emergent 
for  the  earthquake  as  compared  with  the  generadly  impulsive  P  phase  signals  from  the 
explosions.  These  comparisons,  if  universal,  would  make  discrimination  quite  easy  but  a 
casual  inspection  of  Figures  A1-A16  in  the  Appendix  shows  that  this  is  not  imiversally 
true.  For  example  Earthquakes  4  and  5  in  Figures  A4  and  A5  have  relatively  large  P/S 
amplitude  ratios,  more  like  that  of  an  explosion  eilthough  they  still  display  the  emergent 
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P  phases.  For  explosions  in  Figures  A9-A16  the  P/S  amplitudes  are  generally  higher  than 
for  earthquaJces  although  Explosion  1  in  Figure  A 1  is  an  exception.  Explosions  4  and  5 
have  slightly  more  emergent  P  phases  although  the  distinctly  high  single  peak  in  Explosion 
5  indicates  that  the  event  may  be  an  explosion.  The  earthquakes  in  Figures  A1-A8  do  not 
display  such  sharp  initial  pulses. 

The  autoregressive  (based  on  a  third  order  AR  model)  spectra  were  computed  for  the 
P  and  S  phases  for  sill  earthquakes  and  explosions.  This  assumption  is  not  inconsistent 
with  source  models  assuming  a  decay  inversely  proportioned  to  where  v  is  frequency. 
Fitting  a  low  order  AR  spectnim  also  tends  to  de-  emphasize  spurious  peaks  and  values 
due  to  ripple  firing  of  the  mine  explosions.  The  AR  spectra  shown  in  Figures  A1-A8  for 
earthquakes  tend  to  have  stronger  low  frequencies  in  the  S  phases  (0-5  Hz)  and  higher 
frequency  content  in  the  P  phases  (5-8  Hz)  for  earthquakes.  Explosions  generally  tend  to 
have  peaks  at  roughly  the  same  frequency  which  varies  over  the  5-8  Hz  range.  Explosions  1 
and  8  in  Figure  A9  zmd  A16  have  low  to  high  frequency  behavior  consistent  with  that  just 
mentioned  for  earthquakes  so  the  spectral  discriminant  is  not  absolutely  reliable  either. 

To  make  overall  qualitative  assessments  it  is  easiest  to  look  at  the  average  spectra 
of  the  earthquake  and  explosion  groups  separately.  Figvire  2  shows  these  average  spectra 
where  all  traces  have  been  scaled  by  dividing  by  the  maximum  amplitude  of  the  P  phase. 
We  have  plotted  on  a  linear  scale  since  this  emphasizes  spectral  differences  in  the  high  signal 
to  noise  parts  of  the  earthquake  and  explosion  processes.  The  top  left  panel  shows  the 
average  spectra  plotted  on  different  scales  indicated  on  the  left  and  right  ordinates.  This 
display  scales  out  amplitude  differences  and  allow  us  to  see  the  spectral  shape  contrasts. 
For  the  P  phases,  we  notice  a  broader  spectrum  with  both  high  and  low  frequencies;  the 
main  differences  appear  to  be  in  the  low  (0-6  Hz)  band.  We  note  that  in  the  top  right 
panel,  which  plots  the  spectra  on  the  seune  scale,  the  differences  seem  to  be  characterized 
by  a  stronger  low  frequency  (0-6  Hz)  component  for  eeirthquadce  P  thzin  one  sees  in  the 
high  frequency  (6-15  Hz)  bsmd.  The  S  components  shown  in  the  bottom  panels  are  more 
narrow  band  and  relatively  higher  for  eeuihquakes  in  the  interval  (0-3  Hz)  and  higher  for 
explosions  in  the  (3-12  Hz)  band. 

It  seems  clear  that  the  process  of  guessing  spectral  ratios  by  looking  at  separate  or 
average  spectra  will  lead  to  a  niimber  of  possible  discriminants  as  can  be  seen  by  examining 
the  literature.  In  the  next  section,  we  consider  some  amplitude  and  spectral  discriminants 
that  have  been  proposed  in  the  seismic  literature  and  show  that  feat  tires  extracted  from 
the  amplitude  characteristics  do  best  for  the  small  sample  of  earthquakes  and  explosions 
given  above. 
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3.  Discriminants  Based  on  Amplitude  and  Spectral  Features 

For  feature  extraction,  we  consider  a  number  of  classical  measures  related  to  the 
spectrum.  They  are  logarithms  of  (1)  P  and  S  amplitudes,  (2)  P  and  S  mean  square 
error,  (3)  combinations  of  P  and  S  spectra  (0-3  Hz,  3-6  Hz,  6-9  Hz)  zmd  (4)  autoregressive 
coeflficients  of  the  P  and  S  phases  for  a  third  order  model.  The  frequency  ranges  were  not 
exactly  comparable  to  those  used  in  the  literature  (.5-lHz,  2-4  Hz  in  Bennett  and  Murphy, 
1986,  1-2  Hz,  6-8  Hz  in  Taylor  et  al,  1989,  2-5  Hz,  5-10  Hz,  10-20  Hz  in  Pulh  and  Dysart, 
1993,  5-25(5)  Hz  in  Kim  et  al,  1992)  but  were  chosen  by  visually  inspecting  the  separate 
spectra  and  the  average  earthquake  and  explosion  spectra  shown  in  Figure  2.  A  further 
comment  is  that  we  have  avoided  taking  ratios  of  spectra  which  tend  to  assume  a-priori  that 
the  best  discriminator  will  be  the  simple  difference  of  the  form  log(P/S)  =  logP  —  log  5. 
It  is  clear  from  our  results  that  the  log  ratio  is  nearly  the  best  discriminant  and  it  is  also 
reasonable  that  taking  logarithms  improves  the  approximations  to  multivariate  normality. 

The  best  discriminators  of  this  group  were  the  classical  amplitude  and  mean  square 
error  measures  (1)  and  (2);  (1)  is  plotted  in  Figure  3  and  the  scatter  diagram  of  the  mean 
square  error  (2)  hardly  differs  from  this  top  panel.  A  linear  discriminant  analysis  with 
equal  prior  probabilities  tended  to  confirm  the  ratio  procedure.  For  example  the  optimal 
linear  discriminant  functions  for  P  and  5  amplitudes  and  mean  square  errors  were 

d(i)  =  -20.59 log P -I- 15.97 logs -I- 14.32 


and 

d(2)  =  -15.021ogP-f- 13.30 logs  -  25.40 

respectively.  Both  (1)  and  (2)  had  perfect  classification  in  the  test  s3unple  and  classified 
the  first  explosion  as  an  earthquake  in  the  holdout-one  procedure. 

Note  that  the  hold-out  procedure  (see  Lachenbruch  and  Mickey,  1968)  gives  reasonable 
approximations  for  the  misclassification  rates  that  would  obtain  when  classifying  a  new 
observation  not  in  the  training  sample.  The  holdout  procedure  estimates  the  discriminant 
function  for  each  observation  with  that  observation  held  out  of  the  training  sample.  The 
linear  discriminant  function  obtained  is  then  applied  to  the  observation  that  was  held 
out.  The  misclassification  rates  for  all  methods  using  the  original  sample  and  the  hold-out 
procedure  are  shown  in  Table  1. 

The  best  of  the  spectral  group  (3),  also  plotted  in  Figure  3,  focusses  on  the  (0-3  Hz) 
frequency  band  where  differences  were  noted  in  Figure  2.  Of  course,  this  is  bound  to  be 
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closely  related  to  P/S  amplitude  (1)  and  the  mean  square  error  (2)  which  are  both  closely 
related  to  the  low  frequency  power.  It  is  not  sxnrprising  that  it  has  the  same  performance 
as  (1)  and  (2),  leading  to  a  linear  discriminant  function  of  the  form 

d(3)  =  - 13.07 log P  +  12.99 logs  -  10.94. 

Note  that  all  discriminant  functions  are  essentially  of  the  form  —  log  P/S.  This  confirms 
the  intuitive  use  of  measures  of  the  form  log  P/S  which  is  common  in  the  seismology 
literature.  The  -)ther  techniques  in  (3)  (see  Figures  4  and  5)  gave  inferior  performances 
and  Method  (4),  suggested  by  Tj0stheim  (1975)  had  the  worst  performance  with  almost  no 
discrimination  capability.  This  is  surprising  because  the  third-order  AR  spectra  in  Figures 
A1-A16  seem  to  do  a  reasonable  job  of  characterizing  the  spectra  and  because  one  might 
expect  from  standard  source  theory  arginnents  that  such  a  process  would  fit  the  data.  In 
general,  the  third-order  AR  predictions  miss  the  impxilsive  excursions  of  the  explosions.  We 
have  not  investigated  combining  more  than  two  spectral  discriminants  because  of  the  small 
training  samples  (  8  earthquakes  and  8  explosions)  involved  in  the  comparisons.  GlobaJ 
frequency  discriminants  have  been  considered  in  the  literature  for  larger  samples  by  PuUi 
and  Dysart  (1992),  Taylor  et  al  (1989)  and  Kim  et  al  (1992).  Such  global  discriminants 
(see,  for  example,  Kim  et  al,  1992)  can  often  lead  to  linear  combinations  with  both  positive 
and  negative  coefficients.  One  is  hard  put  to  develop  an  intuitive  rationale  for  using  such 
combinations. 

A  comparison  of  performances  on  the  small  test  SEimple  of  Scandinavian  events  is 
given  in  Table  1.  We  include  the  measures  based  on  likelihood  and  information  theoretic 
arguments  given  Section  4  and  the  complexity  approach  given  in  Section  5.  We  note  the 
slight  overall  superiority  of  the  global  optimality  measiures  and  the  general  satisfactory 
performzince  of  the  P/S  amplitude  based  discriminants  considered  in  the  literature.  The 
specific  spectral  ratios  do  less  well  although  it  is  clear  that  the  optimal  tuning  to  the 
spectra  represented  by  the  likelihood  and  information  theoretic  methods  can  do  very  well 
indeed. 
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Table  1:  Misclassifications  (*  =  Holdout) 


Method 

EQ 

EXP 

EQ* 

EXP* 

Amplitude 
logio  P,  logio  S 

0 

0 

0 

1 

Spectral  Discriminants 

MSE  (logio  P,  logio  5) 

0 

0 

0 

1 

•f*(0"3),  logio  *5(0-3) 

0 

0 

0 

1 

logio  P(6-9),  logio  *5(6-9) 

1 

2 

2 

3 

l®Sio  *^(ll'3),  logio  *5(3-6) 

1 

4 

0 

4 

Optimal  Quadratic 

0 

0 

0 

0 

MDI 

0 

0 

0 

0 

ZT 

0 

0 

0 

0 

Complexity 

Bi,d2 

2 

2 

2 

2 

4.  Optimal  Spectral  Discriminants 

For  our  optimal  nonparametric  classification  procedure,  consider  the  classical  approach 
to  discriminating  between  two  stationary  bivariate  Gaussian  processes  {H\  ;  Eaxthquedces 
and  H2  :  Explosions)  with  unequal  matrix  covariance  functions  (spectra).  An  approxi* 
mation  (see,  for  example,  Shumway,  1988)  to  the  optimal  test  statistic  is  related  to  the 
match  between  the  Fourier  spectrum  of  the  series  of  unknown  origin  and  the  spectrum  of 
the  earthquake  or  explosion  process.  Consider  the  likehhood  or  matching  function  under 
Hj,j  =  1,2,  given  by 


dj. 


1 

2 


log/,.(i/fc)  + 


|X(^) 
fi-M  r 


(1) 


where  we  replace  •  by  P  or  5  depending  on  the  phase  to  be  considered,  X.(k)  is  the 
discrete  Fourier  transform  of  the  data  X(.  and  fj.{vk)  denotes  the  spectrum  for  pheise  • 
under  hypothesis  Hj.  The  frequencies  are  of  the  form  i/jt  =  k/T,k  —  0, ...  ,T  —  1.  The 
optimeil  statistic  for  testing  whether  the  sampled  bivariate  series  is  from  Hi  :  EQ  or  from 
H2  :  EXP  is  given  by 


Q  =  {dip  —  d2p)  +  {dis  -  d2s), 


(2) 
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where  we  accept  Hi  (earthquake)  if  Q  >  0  anJ  H2  (explosion)  if  Q  <  0.  Note  that  the  P 
and  S  phases  need  to  be  uncorrelated  processes  for  this  detector  to  be  optimal.  We  have 
computed  cross  spectra  and  coherence  functions  of  the  paired  phases  for  every  event  and 
they  axe  not  significantly  greater  than  zero. 

We  may  also  look  at  information  theoretic  approaches  to  discriminating  between  two 
processes.  It  is  known  (see  Shumway  and  Unger,  1974)  that  the  discrimination  information 
for  two  Gaussian  processes  differing  only  in  the  spectra  is  approximately 


/(/l-,/2.) 


1  ^  r  fii^k) 

2  l/2-(^'*) 


-  log 


f2-M 


(3) 


Generally,  it  is  convenient  to  regard  the  quantity  given  above  as  a  measure  of  the  dis¬ 
crepancy  between  the  two  spectral  densities  /i.(t'fc)  and  /a  (t'fc)  since  /(/i.,/2  )  >  0  with 
equality  if  and  only  if  fi.(vk)  =  f2  (^k)  ^or  all  k.  Kullback  (1978)  has  developed  the 
minimum  discrimination  information  (MDI)  criterion  as  a  means  for  classifying  a  new  ob¬ 
servation  into  Hi  or  H^.  Under  this  principle,  one  compares  the  discrepancy  of  a  spectral 
estimator  computed  from  the  sample  realization  r*.,  say  fr  i^'k)  with  fifi'k)  and  f2  {vk) 
using 


I(h  ,h  -,h  )  =  I(h  ,h  )  -  Hfr  Jy). 


(4) 


where 


nfT.Ji) 


1  fr-ii^k)  /t  (I'fc) 


(5) 


Since  we  want  the  discrepancy  between  the  sample  spectrum  and  the  true  density  to  be 
minimized,  it  is  clear  that  we  should  accept  Hi  when  /(/i  , /2.; /t  )  ^  0  and  accept  H2 
otherwise.  In  terms  of  the  overedl  criterion,  expressed  in  terms  of  both  phases  we  would 
choose  Hi  (earthquake  )  when 


/(/i.,/2.;/Tp)  +  /(/i.,/2.;/rs)  >  0.  (6) 

Note  that  for  fr-i^k)  =  |-X^  (^)Pi  the  periodogram  estimator,  the  above  criterion  reduces 
exactly  to  the  quadratic  criterion  defined  in  Equation  (1).  Zhang  and  Taniguchi  (1992) 
have  shown  the  asymptotic  normality  of  the  MDI  criterion  and  that  the  misclassification 
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errors  converge  to  zero.  They  have  also  shown  that  the  criterion  is  robust  to  departures 
from  normality. 

Zhang  and  Taniguchi  (1993)  have  also  suggested  the  a-entropy 

eM,;M  =  2  E 

0  <  a,  1)  as  an  alternative  and  show  that  it  is  robust  to  both  non-Gaussian  departures  and 
pealc  contamination.  Under  this  suggestion,  we  would  accept  Hi  when  f2  ;fT-)  >  0, 

where 


BaUlJ2\fT)  =  ea(/2.,/r)  “  fx) 


(8) 


M  ] 

m!' 


(9) 


In  terms  of  the  overall  criterion  involving  both  phases,  we  would  accept  Hi  when 


Baifip,  f2P',  frp)  +  Baifis,  f2S',  Its)  >  0.  (10) 

In  order  to  apply  the  discriminant  functions  defined  above,  we  need  to  have  estimators 
for  the  earthqualce  and  explosion  spectra,  say  fi.{v)  and  f2  {v)-  These  can  be  talcen  as 
predefined  values  if  no  training  sample  is  available  or  as  the  averages  of  the  ezirthquake 
and  explosion  spectra  respectively  if  a  training  seunple  is  available.  We  take  the  values 
here  of  the  average  earthquake  emd  explosion  spectra  shown  in  Figure  2.  The  spectra  were 
computed  for  each  series  (no  taper)  over  a  fairly  broad  band  (2  Hz)  and  then  averaged 
separately  for  earthquakes  and  explosions.  Note  that  the  P  and  S  components  were  scaled 
by  dividing  by  the  maximum  of  the  P  component.  For  the  quadratic  and  information  the¬ 
oretic  detectors,  small  values  of  the  theoretical  spectra  can  cause  potential  distortions,  so 
severzil  cutoff  frequencies  in  Equations  (1),  (5)  and  (9)  were  tried;  overall  best  performance 
seemed  to  be  attained  with  a  cutoff  of  about  8  Hz.  In  Figure  1,  we  see  that  real  differences 
in  the  earthquake  emd  explosion  spectra  are  small  after  this  point. 

The  P  and  5  components  of  the  quadratic  detector  are  plotted  in  Figure  6  and  we 
see  that  both  the  test  sample  and  the  holdout  procedures  achieve  perfect  cleissification. 
Note  that  the  performance  of  the  holdout  procedure  emulates  the  performance  that  a 
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discriminant  function  defined  from  a  training  sample  would  have  on  a  new  observation. 
One  computes  the  average  of  the  spectra  holding  out  one  observation  at  a  time  and  then 
using  the  test  statistic  to  classify  the  held  out  observation.  The  events  that  are  somewhat 
"aarginally  classified  since  they  lie  somewhat  close  to  the  decision  fine  in  the  hold-out 
sample  axe  Earthquake  4  (Figure  A4)  and  Explosions  1  cind  8  (Figures  A9  and  A16). 

The  P  and  5  components  of  the  information  theoretic  based  discriminemts  are  shown 
in  Figure  7  for  a  =  .7  eind  we  see  that  the  separation  is  shghtly  better  for  the  robust  ZT 
a- entropy  than  for  the  MDI  detector  whose  performance  is  similar  to  that  of  the  ordinary 
likelihood  discriminant  shown  in  Figure  6.  Experimenting  with  lower  values  of  a  showed 
a  slight  degradation.  It  would  appear  that  the  ZT  a- entropy  discriminant,  which  is  robust 
to  peak  contamination,  is  doing  a  better  job  in  this  case.  Note  also  that  the  distribution  of 
the  ZT  discriminant  is  less  skewed  than  that  of  the  MDI  or  likelihood  detectors  where  the 
explosions  tend  to  be  clustered  in  a  small  region  and  the  earthquakes  tend  to  be  distributed 
over  a  large  dynamic  range.  The  holdout  performance  of  both  detectors,  shown  in  Figure 
8,  is  perfect  and  we  note  that  one  explosion  in  the  MDI  holdout  population  moves  quite 
close  to  the  discriminant  line.  The  marginal  events  might  be  teJcen  8is  Eemthquakes  4  and 
5  although  they  are  quite  far  from  the  decision  line. 

5.  Parametric  Discriminants  Based  on  Complexity 

In  order  to  define  an  optimal  parametric  discriminant,  consider  a  model  for  the  earth¬ 
quake  or  explosion  P  phase  specified  as  a  stationary  autoregressive  series  modiilated  by  a 
time  varying  fimction  sometimes  used  for  eeirthquake  and  explosion  sources.  That  is,  we 
assume  the  observed  P  phase  is  generated  by 


Vi  =  o.t{Q)xt  +  vt  (11) 

where  a((0)  is  some  modulating  function  depending  on  t  and  the  parameter  vector  0  = 
■  The  process  xt  is  an  underlying  signeil  and  its  reflections  and  is  an 
additive  white  noise  process.  Dargabi-Noubary  (1992)  heis  suggested  such  a  model  where 
the  modulating  function  might  be  taken  as 


at{9i,$2)  =  0itexp{—92t)  (12) 

One  may  motivate  such  modulating  functions  by  appealing  to  standard  source  theory 
models  such  as  Harkrider  (1976)  or  Von  Seggern  and  Blandford  (1972)  whose  model  implies 
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a  source  time  function  of  the  form 

at(^i,^2,^3))  =  (^it  +  ^2<^)exp(-^3t). 

The  modulating  function  (12),  suggested  by  Dargahi-Noubzuy  (1992),  has  a  shape  which 
depends  on  the  parameters  6i  and  62.  Small  values  of  61  and  02  should  be  associated 
with  earthquakes  since  they  produce  emergent  modulating  fimctions;  large  values  of  these 
two  parameters  would  characterize  explosions  since  these  laurger  values  will  produce  rather 
impulsive  waveforms.  Some  families  of  typical  modulating  functions  obtained  with  this 
data  axe  shown  in  Figure  9. 

Since  the  modulating  functions  are  rather  smooth,  the  underlying  process  should 
be  modeled  by  a  random  series  with  a  fairly  well  defined  peak  spectrum.  Second-order 
autoregressive  series  are  useful  for  fitting  these  kinds  of  series  and  accordingly,  we  take  the 
modulated  process  Xt  as 

xt  =  +  <f>2Xt-2  +  u;i,  (13) 

where  the  noise  processes  errors  have  variances  <rl  and  for  identifiability,  cr^  should 
be  fixed  at  a  constant  value.  The  second-order  model  implies  a  fall  off  in  frequency  that 
is  inversely  proportional  to  which  is  consistent  with  the  Von  Seggem  Blandford  theory. 

In  order  to  get  an  indication  as  to  how  the  model  defined  in  (12)  and  (13)  might  work, 
we  developed  a  maximum  likelihood  procedmre  for  estimating  the  parameters  01,02,  and 
<j>i,<f>2‘  The  model  is  highly  nonlinear  in  all  parameters  but  we  can  write  the  log  likelihood 
function  of  the  complete  data  eis  in  Shumway  (1988)  and  then  use  the  EM  algorithm.  The 
brisic  procedure  is  to  update  , 4>2  and  cr^  using  the  EM  algorithm  and  to  update  0i  and 
02  by  Newton-Raphson  interations  nested  within  the  EM  algorithm. 

The  above  procedure  is  quite  sensitive  to  start  points  and  the  families  of  modulating 
functions  shown  in  I'igure  9  for  earthquakes  and  explosions  are  not  clear  emergent  and 
impulsive  as  they  should  be.  The  scatter  plots  of  the  parameters  4>2,<t>2  and  0i,02  are 
shown  in  Figure  10  and  agadn  there  is  a  clear  separation  only  relative  to  which  is 
proportional  to  the  P/S  amplitude  ratio.  This  happens  because  the  amplitudes  of  the 
P  is  scaled  by  dividing  by  the  maximum  amplitude  of  the  5  phase  for  that  event.  A 
discriminant  analysis  using  the  parameters  0i  and  02  lead  to  misclassifying  Ezirthquakes  2 
and  8  and  Explosions  4  and  5.  Note  in  Figure  5  that  these  are  the  events  that  one  sould 
expect  to  misclassify  on  the  basis  of  the  estimating  modulating  functions.  Looking  at  the 
original  events,  it  is  reasonable  that  Explosions  4  and  5  would  be  fitted  well  by  emerging 
modulators  but  the  emergent  behavior  of  Earthquakes  2  and  8  would  seem  to  contradict 
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their  fitted  waveforms.  It  is  plausible  that  the  estimation  procedure  could  be  tuned  to  the 
process  by  compsuing  against  envelope  functions  starting  from  a  fixed  time  point  for  the 
maximum  excursion  and  we  are  in  the  process  of  testing  this  method. 

5.  Conclusions  and  Reconrunendations 

We  conclude  that  the  optimal  nonparametric  procedures  based  on  spectral  differences 
discriminate  significantly  better  than  those  based  on  extracting  simple  features  of  the 
process  or  on  fitting  the  amplitude  modulated  model  for  complexity.  Of  course,  these 
results  are  only  for  the  very  small  and  carefully  selected  learning  sample  of  Scandinavian 
earthquakes  zind  explosions  considered  in  this  study.  Hence,  the  data  are  not  sufficient 
to  give  confidence  from  a  discrimination  point  of  view  but  they  eire  adequate  to  indicate 
the  potential  of  the  new  statistical  methods.  It  is  of  potential  interest  also  to  develop  a 
method  for  incorporating  a  third  noise-only  hypothesis  into  the  decision  procedure  in  order 
to  decide  whether  there  is  significant  signal/noise  to  justify  a  discrimination. 

The  advantage  of  the  nonparametric  procedures  based  on  likelihood,  minimum  dis¬ 
crimination  information  and  a- entropy  essentially  relate  to  their  ability  to  tune  against  all 
differences  present  in  the  earthquake  and  explosion  spectra  and  not  to  specific  frequency 
bands  or  phases.  Furthermore,  the  a- entropy  modification  seems  to  be  a  promising  robust 
technique  in  this  case  where  there  can  be  interfering  or  slightly  offset  peaks  in  the  sample 
spectrum  associated  with  the  event  to  be  classified.  In  addition,  for  more  realistic  lairger 
samples  and  more  thaui  one  station  used  per  event,  the  procedure  can  work  even  better. 
The  application  of  the  nonparametric  procedure  to  larger  data  bases  should  provide  some 
standard  baseline  statistics  for  comparison  with  more  esoteric  nonlinear  methods  such  ais 
classification  trees  or  neural  nets  as  applied  to  specific  feature  vectors. 

The  modeling  of  complexity  using  the  waveform  model  proposed  here  may  also  add 
discrimination  capability  for  the  cases  where  the  spectral  matching  fimction  given  by  the 
nonparametric  procedure  is  not  enough.  We  will  continue  to  refine  the  properties  of  the 
waveform  comparison  test  and  apply  it  to  the  small  test  sample  of  Scandinavim  events. 
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Figure  2:  Average  Earthquake  and  Explosion  Spectra.  The  folding  frequency  is  20  Hz  (.5  cycles 
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Figure  3:  Scatter  diagrams  showing  the  separation  between  earthquakes  (open  circles) 
and  explosions  (solid  circles)  for  log  amplitudes  of  P  and  S  (top  panel)  and  log 
spectra  from  0-3  Hz  (bottom  panel). 
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Figure  4;  Scatter  diagrams  showing  the  separation  between  earthquakes  (open  circles) 
amd  explosions  (solid  circles)  for  log  spectra  of  P(0-2  Hz  vs  2-4  Hz)  and  S(0-2 
Hz  vs  2-4  Hz). 
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Figure  5;  Scatter  diagrams  showing  the  separation  between  earthquakes  (open  circles) 
emd  explosions  (solid  circles)  for  log  spectra  of  P(6-9  Hz)  vs  S(6-9  Hz)  and 
S(0-3  Hz  vs  3-6  Hz). 
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Figure  6:  Separation  achieved  for  earthquakes  (open  circles)  and  explosions  (solid  circles) 
using  the  optimal  quadratic  detector.The  top  panel  gives  the  results  using  the 
learning  sample  whereas  the  bottom  shows  the  holdout  scores.  The  separation 
line  Q  =  0  is  shown  in  both  figures. 
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Figure  8;  Separation  achieved  for  earthquakes  (open  circles)  and  explosions  (solid  circles) 
using  holdout'one  versions  of  the  MDI  discriminant  (upper  panel)  and  the  ZT 
discriminant  (lower  panel).  The  sep2iration  line  Q  =  0  is  shown  in  both  figures. 
Note  that  neither  statistic  misclassifies  2uiy  earthquake  or  explosion. 
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Figure  9:  Plots  of  the  estimated  functions  a, (^1,^2)  =  ^itexp(— 52<)  for  the  8  earthquakes 
(Top  Panel)  and  8  explosions  (Bottom  Panel). 
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Figure  10:  Scatter  diagram  comparing  values  of  the  parameters  ^1,^2  corresponding  to 
the  autoregressive  part  of  the  state  space  model  and  and  log  ^2  corresponding 

to  the  modulating  functions  in  Figure  7.  Note  that  most  of  the  separation  is 
with  respect  to  vdues  of  log?i  which  is  proportional  to  log-amplitude. 
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Figure  A2:  Earthquake  2  at  Station  ARAO  on  8/24/91  with  local  magnitude  3.18.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A3;  Earthquake  3  at  Station  NRAO  on  9/23/91  with  local  magnitude  3.15.  P  and 
S  phases  extracted  aire  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A4;  Earthquake  4  at  Station  FlAl  on  1/04/92  with  loca)  magnitude  3.60.  P  and 
S  phases  extracted  are  shown  edong  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A5:  Earthqu2tke  5  at  Station  ARAO  on  2/19/92  with  local  magnitude  3-26.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A7:  Earthquake  7  at  Station  NRAO  on  4/14/92  with  local  magnitude  3.38.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A8:  Earthquake  8  at  Station  NRAO  on  5/18/92  with  local  magnitude  2.74.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectrad 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A9:  Explosion  1  at  Station  ARAO  on  3/23/91  with  local  magnitude  2.85.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  All:  Explosion  3  at  Station  ARAO  on  4/26/91  with  local  magnitude  2.95.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A12:  Explosion  4  at  Station  ARAO  on  8/03/91  with  local  magnitude  2.13.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 

38 


Figure  A13:  Explosion  5  at  Station  ARAO  on  9/05/91  with  local  magnitude  2.32.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A14:  Explosion  6  at  Station  FlAl  on  12/10/91  with  local  magnitude  2.59.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectral 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A15:  Explosion  7  at  Station  ARAO  on  12/29/91  with  local  magnitude  2.96.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectraJ 
estimators.  The  folding  frequency  is  20  Hz. 
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Figure  A16:  Explosion  8  at  Station  NRAO  on  3/25/92  with  local  magnitude  2.94.  P  and 
S  phases  extracted  are  shown  along  with  third-order  autoregressive  spectrsd 
estimators.  The  folding  frequency  is  20  Hz  (.5  cycles  per  point). 
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